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STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR 

DEVELOPMENT 



^ BACKGROUND OF THE INVENTION 

Filter coefficients can be optimized according to many 
^ different criterium. For example stop band attenuation, coding 

Q gain, degree of smoothness can be used as the criterium to 

20 optimize a filter. However, solving for these optimal filter 
coefficients requires solving a time domain problem that is a 
non-linear constrained optimization problem that can be 
difficult to solve. In order to optimize these coefficients 
requires solving the non-linear constrained optimization 
25 equations iteratively to arrive at the optimum solutions. 
However, this iterative process is computationally expensive due 
to the difficulty of solving the equations. In addition, if the 
constraints are not satisfied then the solution is not 
invertible . 

30 Therefore, it would be advantageous to provide a method for 

solving for the optimal parameters for a filter that is not 
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constrained and is more efficient computationally than the time 
domain problem. 



BRIEF SUMMARY OF THE INVENTION 
5 A method for finding optimal filter coefficients for a 

filter given an input data sequence and an objective function is 
disclosed. The method includes selecting a wavelet basis having 
k parameters and minimizing the predetermined objective function 
with respect to the k parameters. The wavelet basis is 
10 reparameterized into k/2 rotation parameters and factorized into 
^ a product of rotation and delay matrices. The k/2 rotation 

O parameters are provided for the rotation matrices and a data 

Q?| transform matrix is computed based on the product of the 

^ rotation and delay matrices. The input data sequence is 

ijj 15 converted into transformed data by applying the data transform 



^ matrix to the input data. The Jacobian of the data transform 

nJ matrix and the input data sequence is determined and multiplied 

^ by the gradient vector with respect to the transformed data of 



the objective function. This product is compared to a 

20 predetermined criterium and if the predetermined criterium is 
not satisfied, a new set of k/2 parameter values are provided 
and the gradient descent is continued until the optimal k/2 
parameters are found. The optimal filter coefficients are then 
calculated based on the optimal k/2 parameters. The wavelet 
25 basis may be selected from a wavelet packet library containing 
orthonormal wavelet packet bases, and in which the selected 
wavelet packet basis minimizes a given cost function, which can 
be an entropy function. 

A method for determining filter coefficients to form a 
30 filter to filter data of length N is disclosed, wherein the 
filter coefficients are optimized to minimize an objective 
function that measures a predetermined quality of the signal 
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data. The method includes the steps of providing the number of 
coefficients, k, in the filter and selecting a wavelet packet 
basis. An objective function is provided and the first set of 

k/2 parameter values are also provided. A data transform matrix 
5 is formed as a function of the selected wavelet packet basis and 
the k/2 parameter values. The data is transformed by 

multiplying the data transform matrix with the data, and the 
value of the objective function is calculated using the 
transformed data. The optimal set of values for the k/2 
10 parameters are then found based on the value of the objective 
^ function. 

p In one embodiment, finding the optimal set of values for 

\l 

^ the k/2 parameters includes first determining if the objective 

O function satisfies a first criteria . If the first criteria is 

^ 15 satisfied, the optimal filter coefficients are calculated using 
the current set of k/2 parameters. If the objective function 
fy does not satisfy the first criteria a gradient steepest descent 

V method is used to modify the k/2 parameters to a local minimum 

P of the objective function. The gradient steepest descent method 

20 is one in which the Jacobian of the data transform matrix is 
calculated and multiplied by the gradient of the value of the 
objective function with respect to the transformed data. If the 
product of the Jacobian and the gradient satisfies a 
predetermined criterium, the iterative process is stopped and 
25 the optimal filter coefficients are calculated based on the 
current set of k/2 parameters. If the product of the Jacobian 
and the gradient does not satisfy the predetermined criterium 
the values of the k/2 parameters are updated based on the value 
of the gradient, and the process is continued in an iterative 
30 manner until the criterium is satisfied and the optimal 
parameters found. 
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other forms, features and aspects of the above-described 
methods and system are described in the detailed description that 
follows • 

5 

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING 
The invention will be more fully understood by reference to 
the following Detailed Description of the Invention in 
conjunction with the Drawing of which: 
10 Figs. lA and IB together form a flowchart depicting one 

method of determining optimal filter parameters according to the 
presently disclosed invention; 

Fig. 2 is a diagram depicting a wavelet packet library of 
orthonormal wavelet packets; 
15 Fig. 3 is a block diagram depicting a Given 's rotation 

block; and 

Fig, 4 is a block diagram depicting a cascade of lattice 
filters suitable for use with the method depicted in Fig. 1 

20 DETAILED DESCRIPTION OF THE INVENTION 

A method for determining the optimal filter coefficients to 
form a filter of length M for filtering data of length N is 
disclosed. The filter is optimized to achieve an optimal value 
for an objective function that measures a predetermined quality 

25 of the signal data. 

Fig. 1 depicts a flow chart for determining the optimal 
filter coefficients^ wherein the number of filter coefficients, 
K, is selected, as depicted in step 102. The value of K can be 
selected based on the filter requirements, the system 

30 requirements, or a priori knowledge of the system response, or a 
combination of one or more these criterium. An objective 
function is provided, as depicted in step 103. The wavelet 
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packet basis^ i.e., the decomposition scheme is selected, as 
depicted in step 104, and a set of K coefficients are selected 
forming a predetermined mother wavelet, as depicted in step 106. 
The K coefficients are reparameterized into k/2 parameters that 

5 correspond to the mother wavelet in the lattice decomposition, 
as depicted in step 108. The data transform matrix is formed 
using the decomposition scheme and the reparameterized k/2 
parameters, as depicted in step 109. The data transform matrix 
is factorized into a product of rotation and delay matrices, as 

10 depicted in step 110. The k/2 parameters are applied to the 
rotation matrices and used to provide the numerical values of 
the data transform matrix, as depicted in step 111. The data 
transform matrix is applied to a data sequence producing 
transformed data, as depicted in step 112. The objective 

15 function of the transformed data is calculated, as depicted in 
step 114, and compared to a first criteria, as depicted in step 
116. In the event that the calculated objective function 
satisfies the first criteria, control passes to step 128 and the 
optimal filter coefficients are calculated using the optimal 

20 parameters. In the event that the first criteria is not 
satisfied control passes to step 118. If the algorithm does not 
converge, a control loop can be used to limit the number of 
iterations and a new set of k/2 parameters can be provided and 
the process restarted, 

25 In step 118, The Jacobian of the f actorized-reparameterized 

second filter basis is calculated. The gradient of the 
objective function of the transformed data with respect to the 
transformed data is calculated, as depicted in step 120. The 
Jacobian and the gradient are then multiplied together to form 

30 the gradient of the objective function with respect to the set 
of k/2 parameters, as depicted in step 122. The gradient of the 
objective function with respect to the set of k/2 parameters is 
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compared to a second criteria;, as depicted in step 124. In the 
event that the gradient satisfies the second criteria, control 
passes to step 128 and the optimal filter coefficients are 
calculated from the first set of filter parameters. In the 
5 event that the gradient does not satisfy the criterium;. control 
passes to step 126. In step 126 the value of the k/2 parameters 
are updated as a function of the value of the gradient value and 
control passes to step 111, wherein the intervening steps are 
repeated with the new set of filter parameters. This process is 
10 continued until either the objective function satisfies the 

5^ first criteria or the gradient satisfies the second criteria. 

O Alternatively, the initial values of the k/2 parameters can 

be determined other than by the time domain methods described 

2 above. For example, after executing steps 102, 103, and 104, 



¥.5 



can be determined randomly and used in the optimization process. 

..... 

Alternatively, a user having a priori information or experience 
Jf can select the values of the k/2 parameters accordingly, 

ill 20 Alternatively, step 105 can also be executed and a mother 
wavelet selected and a set of a priori values for the k/2 
parameters can be supplied. In all of these alternative 
embodiments the optimization of the parameters can continue as 
described above with respect to steps 111-126. Step 128 is 
25 executed in the case where a digital filter architecture is 
used. In the case where the optimization and filtering is to be 
accomplished using lattice filters, there is no need to execute 
step 128 to find the coefficients of the digital filters. 

Filters can be optimized according to different criterium. 
30 The different criteriiom can be described according to a 
predetermined objective function. For example, stop-band 

attenuation, coding gain, and degree of smoothness are three 



15 the values of the k/2 parameters can be determined using other 
techniques. In one embodiment, the values of the k/2 parameters 
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possible criterium. Alternatively ,r the predetermined objective 
function can be based on morphological aspects of the signal to be 
filtered. For instance in medical monitoring applications;, an 
objective function may be used to differentiate between normal 
5 breathing patterns and labored or wheezing breathing patterns, or 
between a normal cardiac rhythm and an abnormal cardiac rhythm. 
The objective function is application specific and developed 
according to a set of promulgated system requirements. 

In general, a discrete wavelet transform includes low-pass 
10 and high-pass filters applied to a signal and the output is 
down-sampled by a factor of two. The high-pass coefficients are 
C5 retained and the high-pass and low-pass filters are applied to 

the low-pass coefficients until the length of the residual 
signal is smaller than the length of the filter. Accordingly, 
15 in general the maximum number of times a filter of length M+1, 
followed or preceded by a decimation by a factor of 2, can be 
nJ applied to a data signal of length N+l, where N+1 is a power of 

two, is given by: 



SI 
01 



PI 



20 Qr^^ floor 



log 



(1) 



It should be noted that the data signal does not have to be a 
power of two for the optimal filter coefficients to be found. 
Although, the algorithm will operate more efficiently and faster 
for signal lengths that are a power of two it is not required. 

25 If the signal length is not a power of two, any of the well 
known techniques to extend the data set may be used. For 
example, zero padding where zeros are added to the end of the 
data signal to increase the signal length without interfering 
with the spectral analysis is an option as is the use of 

30 boundary wavelets to analyze the edges of the data signal. 
Other techniques that are known in signal processing may be used 
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depending on the system requirements, the length of the data 
stream to be analyzed and other design requirements. 

Consider a discrete wavelet transform operating on a signal 
with length that is a power of 2 . As an illustration a filter 
of length 4 is operating on a signal of length 8. Let 
x = [xiOXxO),..jc(7)f be the original signal and c(0), c(l), c(2), and 
c(3) and d(0), d(l), d(2), d(3) be the low pass and high pass 
filter coefficients respectively. After applying the two 
filters to the original data, the result is down-sampled by two: 

[o-(0),cr(l),cT(2),cr(3),^(0),m^(2),^(3)f =C,x , (2 ) 

where 

c(3) c(2) c(l) c(0) 0 0 0 0 ' 

0 0 c(3) c(2) c(l) c(0) 0 0 

0 0 0 0 c(3) c(2) c(l) c(0) 

c(l) c(0) 0 0 0 0 c(3) c(2) 

d(3) d(2) d(l) d(0) 0 0 0 0 ' 

0 0 di3) d{T) d{\) d(0) 0 0 

0 0 0 0 d(3) d(2) d(l) d(0) 

d(l) diO) 0 0 0 0 did) d(2) 

The process is repeated on the just the <r's, the low pass 
coefficients : 

C2[o-(0),o-(l),o-(2),o-(3),J(0),«5(l),<J(2),^(3)f , (4) 

where 

"c(3) c(2) c(l) c(0) 0 0 
c(2) c(0) c(3) c(l) 0 0 
J(3) d(2) rf(0) 0 0 

^ J(l) ^/(O) rf(3) rf(2) 0 0 
^ 0 0 0 0 1 0 
0 0 0 0 0 1 
0 0 0 0 0 0 
0 0 0 0 0 0 



0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


1 



(5) 
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These two steps can be combined and the transformed signal y is 

given by: 

If Ci is orthonorraal, then C2 and C are also orthonormal. 
The conditions for orthonormality are given by: 

c(0)^+c(l)'+c(2)'+c(3)'=l 
c(0)c(2) + c(l)c(3) = 0 

, ^/^\2 _ 1 

(8) 



d{Qf + d{\f + dilf + dOf =1 



rf(0y(2) + t/(iy(3) = 0 

c(0)J(0) + c(l)J(l) + c(2)t/(2) + c(3)J(3) = 0 

c(0)c/(2) + c(iy(3) = 0 (9) 
c(2)J(0) + c(3M(l) = 0 

If the c's satisfy (7), then equations (8) and (9) can be solved 
10 using: 

dik) = {-\fcCi-k),k = Q,...,2, . (10) 
As discussed above, the number of times filters of length M+1 
can be applied to a signal of length N+1, where N+1 is a power 
of two, is given by equation (1) . Accordingly, the matrix C can 
15 be the product of Q orthonormal matrices: 

c — r r c c (11) 

where Q<Qi„ax and Qmax is given in equation (1) . 

The above example is intended for illustrative purposes 
only. In general, the constraints on the c's are given by: 



20 Y.^(n)c{n-2k) = S{k),k = 0,...,{N-l)l2 , (12) 
the constraints on the d's are given by 
^ d{n)d{n - 2k) = d{k\ k = 0,..., (iV - 1) / 2 , (13) 

n=2.k 

and the relations between the c's and the d's are given by: 
XI,, c{n)d{n - 2k) = 5{k), k = 0,...,{N -\)I2 ^ ^ ^ ^ 

cin - 2k)d{n) = 5ik), k = 0,...,(iV - 1) / 2 
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10 



15 



20 



and where 5(k) is the Kronecker delta function. As above, when 
equation (12) is satisfied, equations (13) and (14) may be 
solved using: 

d{k) = (-1)* c{N -k),k-^ 0,..., N . (15) 

This leads to a constrained optimization problem that is 
difficult to solve for some applications- The coefficients 
c{0)r c{l), ... c(N) are reparameterized so that the constraints 
in equation (12) are automatically satisfied. For example, 
equation (7) implies that: 

[c(0) + c(2)f + [c(l) + c(3)Y =U (16) 
which is automatically satisfied by setting 

c(0) + c(2) = cos(^i+^2) 
c(l) + c(3) = sin(^i+^2) 

Where the number of parameters have been reduced by one-half, 
from four to two in the illustrative example. Using the 
trigonometric formulas 
cos(a -h j3) = cos{a)cos(j5) - sin(a) sin(>S) 
sin(a + j3) = sin(ci?) cos(j3) + cos(cir) sm(j3) 
leads to solving for the c's using: 

c(0)-cos(6>j)cos((92) 
c(l) = cos(6'i)sin(6^2) 
c(2) = -sin(^i)sin(^2) 

For a length N filter, the orthonormality conditions imply that: 



(17) 



(18 



(19) 



-l2 



_neven 



+ 



nodd 



= 1, 



(20) 



so that the c's may be reparameterized in terms of trigonometric 
functions given by: 
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(iV-I)/2 AiV-l)/2 A 

2; c(2n) = cos 

«=0 \ n=0 J (21) 

(JV-l)/2 nN-l)/2 \ ' 

2; c(2« + l) = sin 

«-o V n=o y 

As can be seen, the right-hand side contains (N-1) /2 terms, 
but the expansions of the left-hand sides using generalized 
trigonometric addition formula contains 2<''"^'''2 terms. Although 
the reparameterized form of the filter uses less terms, 
distributing the trigonometric monomials to the various filter 
coefficients is difficult and complex. 

Distributing the trigonometric monomials to the various 
filter coefficients is accomplished using lattice factorization 
of filter banks. A polyphase matrix of any two channel 
orthogonal filter bank can be factorized as: 
Hf' =p{e,)Hz)pie,)--K{z)p{0^) (22) 
where p(0)eO(2) and 



A(z) = 



1 0 



(23) 



0 z- 

' cos{0) sin(6')' 
-sinC^*) cos(^) 

This procedure leads to a factorization of a wavelet 
transform. This factorization process is illustrated below for 
a filter having six (6) coefficients operating on a signal of 
length eight (8). Although illustrated for this example, the 
process can be extended for filters having a larger n\amber of 
coefficients operating on signals having a longer length. 

Equation (22) can be rewritten as: 
//C-)=i/<--)A(z)^(^,). (24) 
As an illustrative example, a six coefficient filter using three 
angles will be provided. In particular, the polyphase matrices 
are given by: 
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15 



p p 



(25] 



Hi=H': 



,(3) 



where the orthonormality is assured for jff^ , and the remaining 
polyphase matrices derived from H\ when 



cos(^j ) sin(^, ) 
- sin(^i ) cos(^, ) 



(26) 



(28) 



Thus from equation (24), H^p is given by: 

Hl=H\A{2)p{0,), (27) 
which can be written as: 

" + z-^cf cf) + z-^cf > 1 ^ r cos(^, ) sm(^, )iri 0 11" cos(^, ) sin(6', )" 

+z-Vf ^/fUz-^^/fJ L-sin(^i) cos(^i)lo z-'i-sin(^2) cosC^J, 

Multiplying the right-hand side and solving for the value of the 
various coefficients on the left-hand side can be determined 
using the reparameterized coefficients directly by equating 
powers of z. This leads to: 
c^^^ = cos(6'i ) cos(6'2 ) = c'p cos{02 ) 
cf ^ = cos(6'j )sin(6'2 ) = c^'^ sm(6'2 ) 
cf ' = - sin(^, )sin(^2 ) = -cf > sin(^2 ) 
cf = sin(6'i )cos(6'2 ) = cf > 003(0^ ) 

and similarly for the d's: 
dl""^ =-sm{d, )cos(^2 )= cos(^2) 
^ = -sm(^, )sin(^2 ) = ^^o'^ sin(02 ) 
c/f ^ — cos(^i )sin(^2 ) = -^1'^ sin((92 ) 
Jf> =cos(^,)cos(^2)=ff/'> 005(^2) 

These relationships can be re-written as 



(29) 



(30) 
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C3 






^(2) 










0 






0 



,(1) 



(29) and 



(1) 



0 
0 

sin(^2) 
cos(6'2) 

for the relationship in equation 

0 
0 

sin(6'2) 
003(^2) 

for the relationships in equation 
equation (24), the next polyphase matrix 
function of H^^^ . Multiplying the matrices 
powers of z as above yields: 



(31) 







~ 003(^2) 












0 






0 



d\ 
d^ 



(30) 



H 



(32) 

As noted above in 
can be found as a 
and equating like 







003(^3) 


0 


0 


0 








-sin(^3) 


0 


0 


0 




^^ 




0 


sin(^3) 


003(^3) 


0 








0 


003(1^3) 


-sin(^3) 


0 


c? 






0 


0 


0 


sin(^3) 








0 


0 


0 


003(^3) 





(33) 



Substituting equation (31) into equation (33) yields 











C3 




^2 








.0 . 





003(^3) 


0 


0 


0 






-3in(^3) 


0 


0 


0 


0 


0 


sin((93) 


003(^3) 


0 


-sin(6'2) 


0 


0 


003(^3) 


-sm(^3) 


0 


0 


sm(^2) 


0 


0 


0 


sin((93) 


0 


003(^2). 


0 


0 


0 


003(^3) 







sin(^iy 
cos(^i) 



(34) 



If 5 = ^3^2^ 



then the transpose of B, 



b"^ can be found to be 
Taking the transpose of equation (34) and combining 
equation (34) with a similar matrix equation for solving for the 
d's, yields: 
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(3) (3) 


(3) 

r 

^3 


(3) (3) 
^2 ^1 


(3) 








^ (3) ^ (3) 


^(3) 


^(3) ^(3) 


"o 




COS\^C7j f 




COs(6'3 ) - 




) 0 


0 




0 


0 


0 


0 


• /y-i \ 

sm(^3) 


C0s(^3 ) 


0 


0 


0 


0 


COs(^3 ) 


-sin(6'3) 


0 


0 


0 


0 


0 


0 




sin(^3 ) 


003(^3) 



'cos{0^) -sin(6>2) 0 



0 sin(6^2) cos(^2) 



(35) 

as a function 
Performing the 



of the 
matrix 



Thus, each coefficient can be found 
reparameterized angle coefficients • 
multiplication in equation (35) defines each coefficient as a 
linear combination of two or more trigonometric functions. 
Reformulating equation (4) from above such that the c's and d's 
alternate rows^ the matrix Ci can be written as the product of 
three matrices: 
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o 
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0 0" 
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0 0 
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Where Sk=sin (9^) , Kk=cos {9k) 



10 



(36) 

and all coefficients are assumed to 
be for the third polyphase matrix so that the superscripts have 
been dropped without any loss of generality. The matrix 
equation given in equation (36) can be expressed in the form of: 

c,=R{e,)SRie^)SR(e,) (^v) 

where R(e) are wavelet transform matrices and S are shifting 
matrices. A matrix E is then applied to Ci to separate the high- 
pass coefficients from the low-pass ones. 
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14 10 



15 



20 



25 



The above is for illustrative purposes only^ and signals 
with lengths of 16, 32, 64, or greater powers of two may be 
filtered. As an example, the left-hand side of equation (36) 
can be extended to filter a signal of length thirty-two (32) . 
In this example^ the left-hand side of equation (36) , which is Ci 
in equation (11), becomes a 32x32 matrix in which the 
coefficients are shifted and wrapped around across all thirty- 
two columns, C2 is formed in a block matrix form as follows: 



C 0 
0 ^ 

L ^16x16 



' 16x16 , 



(38) 



where Ci is the left-hand side 8x8 matrix of equation (36) 
extended to a 16x16 matrix by shifting and wrapping the filter 
coefficients as described above. In the illustrative example, 
Qmax is given by equation (1) and is equal to 3. C3 is then 
given in block matrix form as: 



C 0 



0 



16;cl6 



8jc8 

0. 



^ 8;c8 



(39) 



16jc16 ^16;c16. 

where Ci is the left-hand side 8x8 matrix in equation (36) 
Thus, a wavelet transform of data of length 32 can be given by: 



C 0 

^1(8x8) ^8x8 



0 



0 



16x16 



8x8 

0, 



'8x8 



C 0 

^1(16x16) ^16x16 

0 



16x16 



^16x16 J 



[^1(32x32) ] ' 



(40) 



'16x16 16x16. 

This can also be expressed as: 

c^EQRQ{e,)SQRQ{e,\SQRQ{e^y^^^^ (4i) 

As can be observed, the filter is a traditional straight 
wavelet transform- As is known, the straight wavelet transform 
continuously applies a low and a high pass filter to the data, 
down-samples the results by two (2) , retaining the high pass 
results, then applies the high and low pass filters to the low 
pass results only. This is repeated until the filter size is 
too small to filter the remaining data. Alternatively, a subset 
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of wavelet transforms that are less than the maximum number of 
transforms calculated in equation (1) may also be used to form 
the basis. 

Other wavelet bases may be used, however, and the possible 

5 combinations of the placement of the identity matrix and Ci 
matrix within each sub-block can be used. Fig. 2 provides a 
wavelet packet tree corresponding to a wavelet transform of 
input data 202 having a length of eight (8). As discussed 
below, each subspace of the original signal space is filtered 

10 using the same predetermined high and low pass filters. The 
form of the predetermined high and low pass filters is selected 
depending upon the particular wavelet that is being used. The 
input data 202 is filtered by high and low pass filters, G and H 
respectively, into a first subspace having first and second 

15 orthonormal bases 204 and 206 respectively each having a length 
of four (4) . The first and second orthonormal bases 204 and 206 
are filtered by high and low pass filters, G and H respectively, 
forming a second subspace having four orthonormal bases 208, 
210, 212, and 214 each having a length of two (2) . Each of the 

20 four orthonormal bases 208-214 are filtered by a high and low 
pass filter, G and H respectively, forming a third subspace 
having eight (8) orthonormal bases 216, 218, 220, 222, 224, 226, 
228, and 230 each having a length of one (1) . Each possible 
combination of non-overlapping subspaces of the original signal 

25 space forms a unique basis entry in the wavelet packet basis 
library, and in addition, each of the possible combinations is 
also an orthonormal basis. It is possible to identify the 
optimal orthonormal basis within the wavelet packet basis 
library that will provide for an optimal wavelet basis for 

30 representing that particular signal. 

A particularly efficient method of selecting the best basis 
is an entropy-based selection process. In this method, the 
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entropy of a family of wavelet packet bases operating on a 
signal vector is calculated and the best basis is selected as 
the basis providing the minimum entropy. A suitable entropy 
function can be developed in which the sequences of data 
operated on by each basis in the wavelet packet library is 
compared by their respective rate of decay, i.e., the rate at 
which the elements within the sequence becomes negligible when 
arranged in decreasing order. This allows the dimension of the 
signal to be computed as: 



10 = exp 



f \ 

\ n J 

|2(i (|2 



(42; 



where =jx^| j[x|j and -^^Pn^^^i Pn represents the entropy, or 

n 

information cost, of the processed signal. This leads to a 
proposition in which if two sequences {x^] and \x\\ are compared 
so that a function corresponding function and are 

15 monotonically decreasing and if Xo<«<m^" ~22o<«<yK^« then d < d' . 
Although entropy is one measure of concentration or efficiency 
of an expression, other cost functions are also possible 
depending upon the application and the need to discriminate and 
choose between special functions. 

20 The wavelet packet basis selected, which may or may not be 

the optimal basis for the particular signal, is applied to the 
sequence of data forming a transformed sequence of data. The 
particular parameters are optimized to select the optimal set of 
parameters and hence, the optimal set of filter coefficients. 

25 In the preferred embodiment, a gradient descent can be used to 
optimize the parameters. In particular, the gradient of the 
objective function with respect to the parameter set can be 
calculated as: 

V,<Z> = (43) 



-18- 

ATTORNEY DOCKET NO. MEDMT-OOIXX 
WEINGRRTEN, SCHORGIN, 
GAGNEBIN & LEBOVICI LLP 
TEL. (617) 542-2290 
EflX. (617) 451-0313 



where J is the Jacobian matrix for Cx, <j) is the objective 

function, 0 is the set of parameters, and y is the input data 

sequence transformed by the selected basis. The Jacobian matrix 
is of the form: 

Je=[{derC)x {d,,C)x {d^C)x]. (44) 

An explicit form of the derivatives with respect to the new 
parameters can be obtained using equation (41) : 



10 



15 



20 



djC(e,....,e,)=(djCQ)cQ_,...c, +CQ{^CQ_,)..c^ +...CgC2_,...(a,cJ= 

ERq id, )Sq...8jRq (0j ]..S,R, {0, )+ERq [9,)Sq ...djRQ_j0j ]..S,R, {0^ )+ 
....ERq {9, )SQ...d.Rj (Oj \.S,R, {9^ ) = 

ERq(9,)Sq..^J^j +^...S,R,{d^)+ERQie,)SQ.,^Q.i0j +^...S,R,{0^)+ 



ERq{0,)Sq..^^9j +^yS,R,{9^) 



:45; 



where 



8.= — 



j=l,2,...,K and sin' (9) =cos (9) =sin (9+71/2 ) , and 



cos' (6) =sin (6) =cos (9+7r/2) • It is possible to express the above 
derivatives in terms of rotations of the original angles, thus 
making the storage and computation of the derivatives more 
efficient. The derivative of a rotation block is: 



-(a A 




'0 


r 










I 2) 




-1 


0 



m 

so that 

^,RXe^=D,RX9^f 

where Di is a block diagonal matrix, 
becomes : 

d^cie, ,...e^ ) = eRq {e, )Sq ...d^r^ (Oj ]..s,r, {e^ )+ 
eRq {e, )Sq ...De_,i?e.i {e, )..s,R, {e^ )+ 

ERq{9,)Sq...D,R,[9j)..S,R,{9^) 



(46) 



(47; 



Thus, equation (45) 



(48) 
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Accordingly, with the appropriate choice of the matrices Ci, C2, 
. . . f Cq the above expression allows for the computation of the 
gradient with respect to the angles for any basis in the wavelet 
packet library. 

5 Thus, the Jacobian can be computed and stored efficiently 

and separately from the objective function gradient that is 
calculated. Accordingly, the Jacobian is efficiently computed 
for each set of parameters that are used in the gradient descent 
method. For each of the sets of variables that are used, the 

10 Jacobian can be calculated and stored according to equation 
(48) . 

The gradient of the objective function <|> is calculated with 
respect to the transformed data sequence, i.e., the resultant 
data sequence from applying the selected wavelet basis to the 

15 input data sequence. This result is multiplied by the Jacobian 
matrix and the product of the Jacobian matrix and the gradient 
vector forms the gradient vector of the objective function with 
respect to the angle parameters. The gradient vector of the 
objective function with respect to the angle parameters provides 

20 a vector that is the steepest descent from the starting point to 
a local minimum. If the gradient vector is equal to zero, or 
less than a gradient threshold level, a local minimum has been 
reached, or nearly reached, and the objective function is 
calculated for the current angle parameters. If the objective 

25 function is less than a predetermined objective function 
criterium, the current angle parameters are the optimal angle 
parameters and the filter coefficients can be calculated as 
described above. If the gradient vector of the objective 
function with respect to the angle parameters does not indicate 

30 that it is at or near a local minimum, or if the objective 
function exceeds the predetermined objective function threshold, 
then a predetermined adjustment is made to the current set of 
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angle parameters. This process is repeated until an optimal set 
of angle parameters is found. The actual value of the gradient 
threshold levels the predetermined objective function threshold, 
and the predetermined adjustment to the angle parameters are 
dependent on the particular system requirements - 

Fig. 3 depicts a block diagram of a rotation block suitable 
for rotating a pair of input signals according to the rotation: 









«2. 







(49) 



"2, 

where Kg=cos (0) and Se==sin(G). The rotation block depicted in Fig. 

10 3 is considered to have a rotation angle of 0 such that the input 
signals will be rotated through an angle of 0. The rotation block 
300 includes first and second inputs ai and a2. Input ai is 
multiplied by the cos(0i) and sin{0i), and input 0L2 is multiplied 
by the cos(0i) and -sin(0i). The products aicos (0i) and -a2sin(0i) 

15 are added together in summation module 302 to provide output ai^. 
Similarly, the products a2COS (0i) and aisin(0i) are added together 
in summation module 304 to provide output aa^. The rotation block 
depicted in Fig. 3, in combination with delay modules and down- 
sample modules, forms the basic building block of the lattice 

20 filters. When used in a wavelet system, these lattice filters can 
form the high and low pass filters that are used to determine the 
wavelet transform of the input data. 

Fig. 4 depicts one form of a cascade lattice filter that is 
suitable for use with a filter of length four (4) using two 

25 decomposition steps. The architecture and layout of the various 
sub-systems is based on the straight wavelet transform. In this 
formulation the length of the signal can be arbitrary since the 
values of the signal are continuously modified from left to right. 
Fig. 4 is a circuit equivalent of steps 111-118 of the flow 

30 chart depicted in Fig. 1. The lattice filter depicted in Fig. 4 
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is of length four (4) for illustrative purposes only. The filter 
bank depicted therein can be modified to include more parameters, 
i.e., angles, more decomposition steps, i.e., larger Q, and 
different bases from the wavelet packet table by adjusting the 
5 overall architecture in accordance with the system requirements. 
In addition, as the angles are updated during the optimization 
process described with respect to Fig. 1, the corresponding 
rotation angles in the various lattice filters will be updated 
accordingly. 

10 In particular. Fig. 4 depicts a wavelet transform/ Jacobian 

module 400 that receives a sequence of data on line 401 and 
filters the received data with filter module 402. Filter module 
402 is a two stage lattice filter that includes two rotation 
modules as depicted in Fig. 3 that have rotation angles equal to 

15 01 and 02 respectively. Filter module 402 provides a first high 
pass data output and a first low pass data output, 403 and 4 05 
respectively, that represent the high pass and low pass filtered 
input data respectively. The low pass data output 405 is 
subsequently filtered by filter module 404 that includes two 

20 rotation blocks that have rotation angles of 9i and 02 
respectively. The filter module 404 provides a second high pass 
data output and a second low pass data output 407 and 409 
respectively. The first high pass data output 403 and the second 
high pass data output 407 and the second low pass data output 409 

25 together form the transformed data output values y{0), y(l), and 
y (2 ) respectively. 

The Jacobian can also be computed simultaneously with the 
transformed data vector y(0), y(l), and y(2). In particular, the 
data on line 413 is negated by module 414 and provided as a first 

30 input to provided module 412. The data on line 411 is provided as 
the second input to module 412. Module 412 in coit±>ination with 
the portion 402A of filter module 402 form a lattice filter module 
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having first and second rotation blocks having rotation angles of 
Gi and 82 respectively. Module 412 provides a first data output 
415 and a second data output 417, wherein the second data output 
417 is subsequently filtered by lattice filter module 410 that 
5 provides a third and fourth data output 419 and 421. 

The high pass filter data output 403 is provided to lattice 
filter module 408 that has two rotation modules that have rotation 
angles of 9i and 62 respectively. Filter module 408 filters the 
high pass data output 403 and provides a fifth data output 423 and 

10 a sixth data output 425. The data on line 427 is provided to 
module 406^ which in combination with the first portion 404'A of 
lattice filter module 404 also forms a lattice filter having two 
rotation blocks that have rotation angles of 61 and 02 
respectively. Module 406 provides a seventh data output 429 and 

15 an eighth data output 431. 

The Jacobian values are formed as linear combinations of the 
various output data. In particular, J{0) is formed as the sum of 
the low pass data filter output 405 that is negated by module 416 
and the transformed output Y(0) combined in summation module 420. 

20 J(l) is formed as the sum of the fifth data output line 423 and 
the transformed output value Y(2), negated by module 418, combined 
in the summation module 422. J (2) is formed as the sum of the 
sixth data output line 425 and the transformed output value Y(l) 
combined in the summation module 424. J (3) is foimed as the sum 

25 of the first data output line 415 and the transformed output value 
Y(0) combined in the summation module 426. J (4) is formed as the 
sum of the third data output line 419 and the seventh data output 
on line 429 combined in the summation module 428. J (5) is formed 
as the sum of the fourth data output line 421 and the eighth 

30 output data line 431 combined in the summation module 430. the 
data that is provided on the outputs J(0), J (1) , J (5) corresponds 
to the type of decomposition matrix that is used. For example, in 
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the straight wavelet transform depicted in equations (38) -(40) the 
various Y output values correspond to the various blocks of the 
block diagonal matrix. In this example using 64 data input 
signals, J(0) provides 32 values;. J(l) provides 16 values;. J(2) 

5 provides 8 values, J (3) provides 4 values, J (4) provides 2 values, 
and J(l) provides 1 value. As can be seen, the values provided by 
each output corresponds to the position of the identity matrix in 
the straight wavelet transform block matrices. 

It is convenient to require that the mean of the high pass 

10 filter used herein be zero to ensure that no DC component passes 
through the filter without significant attenuation. In the time 
domain formulation a zero mean high pass filter satisfies the zero 
mean filter criteria: 

J^i-iyc{N-n)=0, (50) 

15 where c is the coefficient of the high pass filter as discussed 
above. This condition when translated into the lattice angle 
formulation and combined with the inequalities in equation (21) 
becomes : 



^ . . 2 2 



X(-l)"c(^-«)= £c(2«)-£c(2« + l) = cos -sin 



= 0. (51) 



B=0 «=0 n=0 7=1 J \ 7=1 J 

20 The above is true if: 

where k is equal to N/2, and N is the number of time domain filter 

coefficients. Thus, the search for the optimal filter 

coefficients becomes a constrained optimization problem that is 
25 expressed as: 
argmin 



^^^^.^^/4<^(^i--.^J (57) 



The gradient of objective function can easily be expressed on this 
plane, where the objective function is given by: 
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^ 



(58) 



J 



then the gradient of the objective function provided in equation 
(58) is given by: 



5 As can be seen, constraining the high pass filter to be a zero 
mean filter has the advantage of reducing the complexity of the 
optimization problem by one from k to k-1 degrees of freedom. 

Those of ordinary skill in the art should further appreciate 
that variations to and modification of the above-described methods 

10 for providing optimal filter coefficients may be made without 
departing from the inventive concepts disclosed herein. 
Accordingly,, the invention should be viewed as limited solely by 
the scope spirit of the appended claims. 
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